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The division of the near-field and far-field zones for electromagnetic waves is important for sim¬ 
plifying theoretical calculations and applying far-field results. In this paper, we have studied the 
far-field asymptotic behaviors of dipole radiations in stratified backgrounds and obtained a universal 
empirical expression of near-field/far-field (NFFF) boundary. The boundary is mainly affected by 
lateral waves, which corresponds to branch point contributions in Sommerfeld integrals. In a semis¬ 
pace with a higher refractive index, the NFFF boundary is determined by a dimensional parameter 
and usually larger than the operating wavelength by at least two orders of magnitude. In a semis¬ 
pace with the lowest refractive index in the structure (usually air), the NFFF boundary is about ten 
wavelengths. Moreover, different treatments in the asymptotic method are discussed and numeri¬ 
cally compared. An equivalence between the field expressions obtained from the asymptotic method 
and those from reciprocal theorem is demonstrated. Our determination of the NFFF boundary will 
be useful in the fields such as antenna design, remote sensing, and underwater communication. 


1. INTRODUCTION 


In many cases, an electromagnetic field shows near-field or far-field behavior when it is observed in a region near 
or far from the field source. The division of the near field and far field is not only helpful for simplifying theoretical 
calculation in different regions, but also easy to highlight the features of the field in respect regions so as to provide 
physical interpretations of the field behaviors. An important question of how to distinguish the near and far fields 
naturalfy arose. Usually, the field can be expanded by negative powers of the distance between the field point and 
source [ll-l3[. The far field is dominated by the lowest order term, and the near field by the higher-order terms. The 
former arises from the requirement of energy conservation and is call a “far field approximation” or “leading-order 
approximation (LOA)”. 

For dipole radiations in a vacuum background, the near-filed/far-field (NFFF) boundary is of the order of light 
wavelength Aq, L ~ Aq. As the source becomes larger, this dimension should be corrected as L ^ D'^/Xq by the 
diffraction theory, where D represents slit widths or antenna lengths. However, what we frequently encounter in 
practical situations are stratified backgrounds, rather than the vacuum one. Stratified backgrounds are of complex 
configurations, in which the NFFF boundary has not been properly addressed yet. To reveal the NFFF boundary in 
such backgrounds is not only a complement in theory, but also very important for many practical applications, such 
as remote sensing, antenna design, NFFF transformation and so on [lj-l9|. In these applications, the distances between 
sources and observation points are usually much larger than operating wavelengths, so that the sources can be treated 
as dipoles. Here we investigate the far-field asymptotic behaviors of dipole radiations in stratified structures, and 
address a universal NFFF boundary. 

Calculation of the dipole radiations in stratified structures is mathematically equivalent to dealing with Sommerfeld 
integrals (Sis). In general, numerical evaluations of the Sis are not easy since the integrals have an oscillatory feature 
and possess singularities along or near the integration paths. Several numerical methods have been proposed to 
calculate the Sis so far. Paulus et al. [l^ developed a direct numerical integration method by appropriately choosing 
the integration path that avoided all possible singularities in the complex plane. This method has been considered 
as a standard test for the accuracies of other methods. The discrete complex image method [lll - [l^ and the rational 
function fitting method [3 [13 tried to expand the integrands by simple functions to obtain closed-form solutions. On 
the other hand, under the far-field approximation, closed-form expressions could be obtained from Sis by asymptotic 
methods [3-[6l. [l7l - [3lj| or reciprocal theorem 32-^. Compared with the numerical ones, the asymptotic methods are 
of advantages of simplicity, easy programming, and clear physical interpretations. 

The stationary phase method and the steepest descent method are two commonly used asymptotic methods [2“ 
Sj [3 - l3| . Both of them approximate the value of a rapid-oscillating integral by the contributions around stationary 
points (SPs) or saddle points, and provide the same LOA results for the Sis (3- In this paper, we adopt a simplified 
version [20|, [2l| of the stationary phase method to perform the asymptotic analysis. This simplified version has 
been successfull y appli ed in the research of underwater communication 22-24|, antenna design 3 23; NFFF 
transformation 1271 l28l . 


Theoretically, the LOA can acquired precise results from Sis when the observation point moves to infinity [18|, [21 1. 
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But for practical usage, an empirical distance is needed to justify its applicability, which is the NFFF boundary. 
Although an error analysis by estimation of higher-order contributions may reveal the boundary, the analysis becomes 
very difficult in stratified structures. Here we compare the asymptotic results with the accurate numerical ones in 
different stratified structures. In this way a universal empirical boundary is obtained. It is found that the NFFF 
boundary is mainly affected by lateral waves, which correspond to the branch point contributions in Sis [^-Q. Besides, 
the boundary is sensitive to the structure configurations and the location of the dipole. Moreover, two kinds of 
treatments in the asymptotic method are carried out and their numerical results are compared. The equivalence 
between the field expressions obtained from the asymptotic method and reciprocal theorem is also demonstrated. 

The paper is arranged as follow: Section 2 sets our model of the stratified structures and presents the LOA 
formalism. Section 3 and 4 discuss the far-held distributions and the NFFF boundary in bilayer and trilayer structures, 
respectively. Then the results are generalized to multilayer structures and proved to be universally applicable in Sec. 
5. Finally, Sec. 6 presents our concluding remarks. 


2. MODELING THE STRATIFIED STRUCTURE 
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Figure 1: Schematic of a stratified structure. The cylindrical coordinates are set and the origin is always on the VED. 


The model studied in this paper consists of a lossless layered structure with the stratihcation along the z direction. 
It contains Q layers and Q — 1 interfaces from top to bottom, as depicted in Fig. 1. Since we concentrate on the far 
held, the two semispaces, i. e., the top and bottom layers, are focused on and the intermediate regions are ignored 
when their details are not needed. The lower surface of the top layer is at z = di and the upper surface of the bottom 
layer is at z = dg-i, respectively. A vertical electric dipole (VED) is located on the origin inside the structure. For 
convenience, we consider the z component of the electric helds generated by the VED. Other held components 
and other dipole orientations will be discussed in Sec. 5. In Fig. I, the cylindrical coordinates (p, ip, z) are set. If the 
VED is either in the top layer or in the intermediate region, the held in the top and bottom layers is expressed as 


Ez{p,z) 


[P.F.]S,, 




y_oc 


/ oo ^3 

-OO ^ 

^qz 


z > di- 


z < dg-i- 


( 1 ) 


The quantities in this expression are as follows. represents the Hankel function of the hrst kind; jz is the 

dipole moment in the z direction; Ci and Cq are respectively the scattering coefficients in the top and bottom layers; 
kq = uj^EqPqjc = -|- is the wave vector in the gth layer; w is the angular frequency; c is light speed in 

vacuum; p,Q and pq are respectively the permeability in vacuum and in the gth layer; P.F. standing for primary held, 
is the Ez generated by the VED in a homogeneous background; 6iq is Kronecker delta. In Eq. (I), the VED is not 
allowed to appear in the bottom layer. However, if the VED is in the bottom layer, one may use the reversed geometry 
instead. 
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3. BILAYER STRUCTURES 


This section discusses the far-held asymptotic behaviors of dipole radiations and the NFFF boundary in bilayer 


structures. 


Reflection Transmission 



Figure 2: Physical meanings of SPs in a bilayer structure. Method A: blue lines; Method B: yellow lines. For reflections, 
Method A and B show the same picture that the reflected field is equal to a field generated by an image point of VED mirrored 
by the interface. For transmissions, Method A describes a refraction that satisfies Snell’s law, while Method B sets the start 
point on the interface right below the VED. 


In this case, the intermediated region in Fig. 1 is removed, and the interfaces at z = di and dg-i merge into one, 
as shown by Fig. 2. The top and bottom layers are also denoted as layer 1 and 2, respectively. 


3.1. The expressions of the Ez in the far-field zones 


In the case of bilayer structures shown in Fig. 2 (di < 0), Ci and Cq in Eq. (1) are expressed as 


f Cl = 


kiz/si - k2z/£2 
klz/si + kizl£2 
“^kiz/ei 

klz/Sl + k2zl£2 


( 2 ) 


where £q is the relative permittivity of the qth layer. In this paper, and ty denote the Fresnel scattering coefficients 
of a single interface when light incidents from layers i to j. Here we carry out a detailed derivation on the transmission 
field with two asymptotic treatments. In this process, the similarities and discrepancies between these two treatments 
are disclosed. 

Substitution of Eq. (2) into (1) and application of the asymptotic form of Hankel function give the expression of 
Ez in layer 2: 


Ez{z < di) 


/ OO LO 

-OO kiz 


( 3 ) 


The SP is obtained by taking derivative of the phase as follows: 

p+-^di +-^(z - di) = 0. (4a) 

klzs k2zs 

This way determining the SP is called Method A. Its physical meaning can be explained by ray theory. On the other 
hand, if the factor exp(—in Eq. (3) is a slowly varying one, it can be taken out of the integral, and the 
remaining phase gives the SP as follow: 


P+-^{z-di) =0. 
f^2zs 


m 


This way is called Method B [2^, [27|. In Eqs. (4a) and (4b), the subscript “s” stands for stationary points. Here we 
emphasize the features of these two methods. Method A treats exp(—i/ci^di) as an oscillatory factor. As a result, 
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the term kpgdi/kizs is included in Eq. (4a) in determining the SP. It describes the wave paths in both layers and 
gives the refraction picture shown by the blue lines in Fig. 2. By contrast, Method B treats exp(—as a slowly 
varying one so that the term kpsdi/kizs is excluded in Eq. (4b) in determining the SP. This description is merely 
applicable to the wave path in layer 2, as shown by the yellow lines in Fig. 2. 

Under these points of view, Method A gives the electric field as 


Ez{z < di) = - 


87rfc| 


. (i. r hiL 


H^^\kpp)e-^’^^‘ 


{z — d\)—ikizd\ 


(5a) 


while Method B gives 


Ez{z < di) = 


87rfc| 


r^lzS 



dkp^H^^\kpp)e-^^^‘\^-’^^\. 

f^2z 


{5b) 


The simplification in Method B is that the integral in Eq. (5b) can be expressed as spherical waves using the 

Sommerfeld identity @,113 


„ikr 


dkp^H‘'^\kpp)e-^’^^\^\. 

, ^Z 


By contrast, Eq. (5a) cannot be further reduced in mathematics because its exponential term contains two different 
wave vectors. We have to carry out the integral according to Eq. (4a). The corresponding ray interpretation is 
demonstrated by blue lines in Fig. 2. 

After integrations, the asymptotic expressions of E^ in Method A and B are respectively given by 


E:^{z <di) = - 


^k2zs ,^2 

Sirk. 
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(6a) 
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ps 
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( 66 ) 


The expressions in the square brackets in Eqs. (6a) and (6b) are the resultants of the integrals in Eqs. (5a) and 
(5b), and those in the curly brackets represent the propagations of E^ along the blue and yellow lines shown in Fig. 
2, respectively. Thus it is easy to see that Eq. (6a) considers E^ as spherical waves in both layers and describes a 
refraction process at the interface. The factor exp[i(A:i — fe 2 )Ari] in Eq. (6a) is the phase correction for the spherical 
wave propagating in layer 1, as can be seen by comparison of Eqs. (6a) and (6b). The term ti 2 {kps)k 2 zs/kizs arises 
from the boundary condition. 

The physical explanation of Eq. (6b) is similar. The factor exp(—i/ci^di) reflects the influence of the wave in layer 
1. In layer 1, is a plane wave propagating along the z direction, and after crossing the interface at the point right 
below the VED, it turns to the spherical wave and propagates to the observation points. 

Method B was given by Refs. [23 271. Comparatively, Method A provides a more accurate description and 


asymptotic expression for the transmission than Method B. That is why we suggest Method A in this paper. However, 
the discrepancies of these two methods tend to be trivial as di/r goes to zero. Another treatment of the transmission 


field |22h24| employed Eqs. (4a) and (6b). In other words, it considered exp{—ikizdi) first as an oscillatory term to 


calculate the SP, and then as a slowly varying one to take out the integral. The accuracy of this method is between 
Method A and B. We do not discuss it in the following. 























5 


Now let us show the equivalence between the asymptotic method, Method B, and the reciprocal theorem. Eq. (6b) 
can be reformed as 


E^{z < di) 




( 47rfc| r2 j 


(7) 


because 


7:<^ 


£2 ‘^k2zsl^2 

^Izs /ei + ^2zs /e2 


— t2l{kps). 
El 


Since t 2 i{kps) is the transmission coefficient when light incidents from layers 2 to 1, Eq. (7) describes the reversal 
propagation picture: the is now generated by a VED at far-held zone and propagates back alone the reversal 
direction of the yellow line in Fig. 2. After crossing the interface, the spherical wave turns to a plane wave and 
propagates to the origin. That the spherical wave emitted by a test source turns to a plane wave after passing the 
interface is the treatment of the reciprocal theorem. In Ref. [s^, the turning point was chosen such that it was at 
the line along the direction of the VED. Therefore, the result there was the same as Eq. (6b) here. The conclusion is 
that Method B is equivalent to the reciprocal theorem. 

When the observation points and the source are in the same layer, the reflection process described by Methods A 
and B are identical , as shown in Fig. 2. The asymptotic expression of E^ in layer 1 is expressed as 



[P.F.]5iq P ri2{kps) 
+ {z- 2 di )2 


iuJUoUij^ 2 

ri j ■ 


( 8 ) 


3.2. The features of the fields 

With the formulas derived above, we are able to do numerical computation. The numerical results of the asymptotic 
methods are compared to the precise numerical ones 0. This enables us to acquire the NEFF boundary. As an 
example, we take Ag = 1.5/im and /i^ = 1. As can be seen later, our conclusions are valid for other wavelengths and 
permeabilities. Since is of an axial symmetry, the angular direction ip is not considered. Instead, the observation 
points are labeled by the spherical coordinates {r,9,ip). Figure 3 shows the value of |E^| as a function of angle 9, at 
different distance r. When 0 < 0 < tt/2, the observation point is in layer I, while when tt/2 < 9 < tt, the observation 
point is in layer 2. In Fig. 3, the dielectric constants in layers 1 and 2 are set to be I and 2.25, respectively. So, the 
upper and lower layers are called lower and higher refraction index (RI) layers, respectively. 



Figure 3: \Ez\ as a function of angle 9 for different distance r. Structural schematics are illustrated in the insets. The 
critical angle is 9c '~(180°-42°)=138°. (a) (|di|, r)=(10,1000)Ao. (b) (|di|, r)=(10,10000)Ao. (c) (|di|, r) = (0.1,10)Ao. (d) 
(|di|,r)=(0.1,100)Ao. 


In Fig. 3, the precise numerical results 0 are plotted by red lines for comparison. For reflections, the oscillation 
curves in Figs. 3(a) and 3(b) reflect the interference between the primary field and reflected field, while this pattern 
is not shown in Figs. 3(c) and 3(d) because the distance between the VED and its image is too close. The asymptotic 
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Ao [23|, which agrees with the 


results fit the numerical ones very well when r is of the order of IOAq, or smaller r 
NFFF boundary in a vacuum background. 

When transmitting from the lower RI layer to the higher RI one, there exists a critical angle Oc above and below 
which are allowed and forbidden regions, respectively [36|, l37[. The boundaries of the forbidden regions are indicated 
by the transmission peaks in Fig. 3. Please note that 0^ is related to di value, so that the transmission peaks in 
the four panels of Fig. 3 have slight differences. In Figs. 3(c) and 3(d)^he curves in the forbidden regions show 
ripples which is an interference effect caused by the boundary continuity |35| . Since the ripple is a near field effect, it 
does show up in Figs. 3(a) and 3(b) where the observation distance is far. When the distance is less, the ripple will 
become evident. Moreover, when the real part of kiz vanishes, exp(—z/cizdi) should be treated as a slowly varying 
term, which means that Method A can be replaced by Method B near 9c and within the forbidden region. 

In the allowed region, although Method A provides more accurate results, both of them coincide with the numerical 
results. However, within the forbidden region and around 0c, there is a big difference between our results and the 
numerical ones. When the VED is far away from the interface, as shown in Fig. 3(a), the asymptotic result decays 
rapidly around 0c due to the effect of exp(—ifcizdi), while the red line shows a slow attenuation. At first glance, this 
seems strange since the numerical results also contain exp(—ifei^di). This can be answered by Eqs. (3) and (6b). 
The numerical method treats exp(—in Eq. (3) as a whole factor, which means that when both of the 
numerator and denominator tend to 0, kiz slows down the decay of exp(—i/cizdi). While the asymptotic methods 
treat l/fciz as a constant, implying that only the decay of exp(—i/cizdi) is considered. In Sis, fciz = 0 corresponds to 
the fei branch point, and at this point the integration gives a lateral wave Lateral waves are the surface waves 

that mainly exist in the higher RI semispace and exponentially decay in the lower RI one. This explanation is also in 
line with the description of the steepest descent method. Near 0c, the steepest descent path approaches the ki branch 
point and its contribution to the integral becomes inegligible. In the forbidden region, the branch point is enclosed 
within the deformed integration path, and causes the lateral waves to appear in the forbidden region. Since the lateral 
waves represent a higher-order attenuation, it is not covered by the LOA. This indicates that the differences between 
the asymptotic results and numerical results come from the branch point contributions. 



Figure 4: Results of the power series fitting the amplitude of the lateral waves at angles 120° and 135° in Fig. 3. Fitting 
functioin is \Ez\=ar~'^^ where a (not show here) and n are fitting coefficients. The fitting range is r=4000 lOOOOAo [3^ . 


The branch point contributions can also be treated by asymptotic methods, e. g., branch cut integrals based on the 


steepest descent method 
as a slowly varying term [( 


or uniform asymptotic expansions [6(. However, these methods regard exp(—ffcizdi) 
Hlai, is resorted generalize 

such asymptotic expressions to the cases of multilayer structures [^ . The exploration of asymptotic expressions for 
branch cut integrals is beyond the scope of this paper, so here we numerically calculate the decay behaviors of the 
lateral waves at angles 120° and 135° in Fig. 3, and show them in Fig. 4. It is seen from Fig. 4 that when |di| > 2Ao, 
the lateral waves decay as r“^. This coincides with the case shown in Figs. 3(a) and 3(b) where r decreases by one 
order of magnitude while the field decreases by two orders. When Aq < |di| < 2 Aq, the way of the decay of the lateral 
waves turns to so that in Figs. 3(c) and 3(d) the numerical results and the asymptotic results have similar 
distributions in the forbidden region. 


3.3. Determining the NFFF boundary 

Having the knowledge of the decay behaviors of the lateral waves, we discuss the NFFF boundary in the bottom 
layer. As can be seen from Fig. 3, the boundary may be chosen as L ^ 1000|di|, which satisfies the evaluation of the 
order of magnitude indicated by Fig. 4. However, this distance should be modified because of two reasons. One is that 
for a very small |di|, e. g., |(ii| = O.OOlAo, A ~ Aq may not eliminate the influence of the lateral waves (numerically 
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verified), which means that a lower boundary should be better. Considering Fig. 3(d), the boundary may be written 
as L ^ 1000xMax{|di|, O.lAo}, where Max{, } means “choose the larger one”. The other is that compared to the one 
of a homogeneous background, the NFFF boundary here has a looser relationship with Aq. Although the two |c?i| 
values in Fig. 3 differ by 100 times, the NFFF boundary only depends on |di|. Moreover, as implied by Eqs. (3), 
(6a) and (6b), the field distribution has a scale invariance over Aq in a non-dispersive media, indicating that the field 
behaviors shown in Fig. 3 stand also for other light wavelengths. So, what influences does Aq have on the NFFF 
boundary? Except for the lower boundary, it mainly decides the location of the forbidden region and introduces 
dispersion. The dispersion considered, the NEFF boundary is suggested to be 

77 

L - 1000^^ X Max{D, O.IAq}, (9) 

^min 

where Umax a-nd respectively represent the largest and smallest RI in the structure, Uq = (eg/Tg)^/^. Equation 
(9) is the most significant conclusion of this paper, and it is valid for any stratified configuration, as will be illustrated 
below. Here, a dimensional parameter D is introduced. In the present bilayer structure, D equals to |(ii|. Eor other 
configurations with more layers, D needs to be modified, as will be shown below. 

In the following, we test Eq. (9) with other configurations. To avoid repetition, the subsequent discussion will not 
compare the results of different observation distances, but only focus on those at r = lOOAo, a distance less than that 
given by Eq. (9). 



Figure 5: \Ez\ as a function of angle 0 at distance r=100Ao. Structural schematics are illustrated in the insets. The critical 
angle is 9c ~ 42°. (a) |di| = 10Ao. (b) |di|=0.1Ao. The coordinate between 30° ~ 96° is stretched to clearly show the oscillation 
details. 

Now we set the dielectric constants in layers I and 2 to be 2.25 and I, respectively. So, in this case the VED is in 
a higher RI layer. as a function of angle 9 for two |di| values is plotted in Fig. 5. In the range 0 < 0 < 7r/2, i. 
e., in layer I, the oscillations shown in Figs. 5(a) and 5(b) reflect the interference between the primary and reflected 
fields as have been mentioned in analyzing Fig. 3. Since the lateral waves favor to appear in the higher RI semispace, 
in the present case, they mainly influence the reflection. 

In Fig. 5, the coordinates between 30° ~ 96° are stretched to illustrate the oscillations. Around 0c, the asymptotic 
results and the numerical results have slight differences in the range 45° ~ 50°, so that Eq. (9) is also applicable here. 
As for the transmission, since now light is from the layer with lower RI to that with higher RI, no forbidden regions 
appear in this configuration. Light refracts and the asymptotic results fit the numerical results quite well. 

According to the discussions above, it is known that the NFFF boundary is mainly affected by the lateral waves. 
In the higher RI semispace, the boundary is given by Eq. (9) and much larger than the operating wavelength; while 
in the lower RI semispace, the boundary is about ten wavelengths. Please notice that it is not appropriate to define 
the NFFF boundary by distinguishing the allowed and forbidden region, because there is no forbidden region in the 
configurations shown in Fig. 5. 


4. TRILAYER STRUCTURES 

Now we investigate the case of trilayer structure. The three layers are called the top, middle and bottom layers, or 
layers 1, 2 and 3, as illustrated in Fig, 6. Two cases are involved where the location of the VED is in the top and 
middle layers, respectively. 
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Before presenting numerical results, three illustrations ought to be addressed. Firstly, different from bilayer struc¬ 
tures, there will occur in a trilayer structure three kinds of surface modes guided modes, proper complex 

modes and improper complex modes, which correspond to poles in Sis. The first two belong to normal modes that 
build up the point spectrum in eigenvalue problems. While the last one is commonly referred as leaky modes which 
violates the radiation condition and only shows up in a certain range of angle in space [s^, . The guided modes 

occur in the middle layer of a mode guiding configuration where the middle layer possesses the largest RI. The proper 
complex modes cannot propagate far from the source due to their decaying nature. Therefore, these two kinds of 
modes will not be touched in the following discussions. The leaky modes affect the far-held distribution and will be 
considered. Secondly, because there are two interfaces in this structure, light multi-reflects in the middle layer and 
may cause more oscillations in the far-held distribution. We will not explain all the formations of these oscillations, 
but focus on the ones that may affect the NFFF boundary. Moreover, the differences between Method A and B are 
not obvious in the previous section because dxjr are small, but will be magnihed by the multi-rehection in the middle 
layer. Thirdly, there are six possible conhgurations of dielectric constant distribution for a trilayer structure. We 
merely study two of them: [e : 1 — 2.25 — 12] and [e : 2.25 — 12 — 1], which are often encountered in experiments. Our 
conclusions are easily extended to other conhgurations. 


4.1. The VED is in the top layer 


Reflection Transmission 



Figure 6: Physical meanings of SPs when a VED is located in the top layer. Method A: blue lines; Method B: yellow lines. 
Concerning reflections, Method A shows that the multi-reflection process in the middle layer contributes to the reflection in 
layer 1, while Method B only considers the direct reflection from the interface at d\. For transmissions, Method A shows that 
the multi-reflection process contributes to the transmission in layer 3, while Method B sets the start point on the interface at 
d 2 right below the VED. In the figure, only the modes m=0 and 1 in the multi-reflection are depicted. 

When the VED is in the top layer, as shown in Fig. 6 (^2 < d\ < 0), the scattering coefficients Ci and Cq are 
expressed as 




^—ikizdi 


di 


( 10 ) 


where the expressions of the coefficients rij and tij can be found in Eq. (2). After expanding Eq. (10) with the 
geometrical optics series and following the steps in Sec. 3, the SPs of Method A for reflections can be obtained by. 


P- 


l(™) 

^ps 

l(™) 




2di) 


Am) 

' 2m{di 
Am) '■ 


d2) = 0, 


( 11 ) 


where m G [0,oo) is the order of the geometrical optics series. The SP of Method B is obtained by letting m = 0 in 
Eq. (11). In the same way, the SPs of Method A for transmissions meet 


P + 


Am) 

^ps 

l("^) 

^Izs 


di 


^(2-+l)(rfi 




(m) 


Am) 

Azs 




d2) = 0 . 


( 12 ) 
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Since Method B only consider the phase variations in the observation layer, its transmission SP has the same form as 
Eq. (4b), with di and k 2 zs being replaced by d 2 and respectively. The physical meanings of the SPs are depicted 
in Fig. 6. Method A describes a multi-reflection process, so that its expressions for reflections and transmissions 
are in the form of series expansions. Each term of the expressions is similar to Eq. (6a), where Ar 2 now represents 
the multi-reflection path in the middle layer and it is also needed for describing the path in layer 3. As for Method 
B, the asymptotic expressions are similar to Eqs. (6b) and (8). Moreover, the equivalence between the expressions of 
Method B and reciprocal theorem can also be proven, in a similar way to that done in Eq. (7). 



Figure 7: \Ez\ as a function of angle 9 for four geometry configurations with [e : 1—2.25 — 12]. Two critical angles are 6 ci ~(180° — 
17°)=163° and 0c2 ~(180° — 26°)=154°. Structural schematics are illustrated in the insets, (a) (jdij, |d 2 |, r) = (0.1,1.1, 100)Ao. 
(b) (jdij, |d2|,r)=(0.1,10.1,100)Ao. (c) (jdij, |(i 2 |, r)=(l.l, 10.1,100)Ao. (d) (jdij, |d 2 |, r)=(0.1,10.1, 35000)Ao. The shaded areas 
between 6 ci and 9 c2 are the second total internal reflection regions. 

Figure 7 plots the distributions in the case of [e : 1 — 2.25 — 12] configuration. Since layers 2 and 3 have higher 
RI, there are two critical angles: 9ci ^ — asin{l/y/l^) and ~ tt — asm(l/V2.25). 9ci determines the forbidden 

region. The range between 9ci and 9 c 2 is shaded in Fig. 7 and is termed as the second total internal reflection region. 

In Fig. 7(a), the thickness of the middle layer is Aq. It is seen that the differences between the asymptotic results 
and numerical results emerge in the transmission field. These differences are caused by the lateral waves as well. The 
dimensional parameter D in Eq. (9) now should equal to the distance between VED and the lower interface at z = d 2 . 
The calculated NEFF boundary is approximately 35 times of the observation distance, which is enough to eliminate 
the influence of lateral waves (numerical verified, not show here). 

In Fig. 7(b), the middle layer is thicker: its thickness is IOAq. This figure shows two significant features. One is that 
Method A and B have noticeable differences in both of the reflections and transmissions. Because Method A takes into 
account the multi-reflection, its results are closer to the numerical ones. The other is that in the second total internal 
reflection region, the asymptotic curve shows rapid oscillation [The coefficients Ci and C 3 in Eq. (10) are of the form 
of multi-reflection, and in this sense, there is also multi-reflection in Method Bj. The formation of the oscillation can 
be explained as following: although the light decays in the top layer, part of the energy can be transferred into the 
middle layer through evanescent waves and become propagation ones. Then in the middle layer it scatters at the 
lower boundary, and has a total internal reflection at upper boundary, which further forms a multi-reflection to yield 
the oscillation. When the VED moves away from the interface, the oscillation will die away. This is verified by Fig. 
7(c) where jdij changes from O.IAq to Aq. It is seen from Fig. 7(c) that the oscillations disappear and the far-held 
distribution decreases. It is worth mentioning that the description here is the picture of leaky waves. A comparatively 
thicker middle layer allows many leaky modes to exist, resulting in the rapid oscillation in the far-held distribution. 
As a test of the NEFF boundary. Fig. 7(d) plots the held distribution of the same conhguration in Fig. 7(b) but with 
the observation point being at the NFFF boundary, i. e., with the observation distance computed by Eq. (9). The 
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results fit each other very well. 



Figure 8: \Ez\ as a function of angle 9 at distance r=100Ao for two geometry configurations with [e : 2.25 — 12 — 1]. The two 
critical angles are 6ci ^ 26° and 9c2 ~ 42°. Structural schematics are illustrated in the insets. |di|=0.1Ao. (a) |d 2 |=l.lAo. (b) 
|d 2 | = 10.1Ao. The fields in the middle layer are not presented. 


Fig. 8 shows the \Ez\ distributions in [e : 2.25 — 12 — 1] configuration. Since the middle layer possesses the largest 
RI, the structure supports the guided modes. The peaks around 90° of the numerical results are the decays of the 
guide modes in layers 1 and 3. Moreover, there is no forbidden region in the structure. Consequently, the asymptotic 
result and the numerical results fit very well. For reflections, these two kinds of results show some differences around 
9 1-2 due to the effect of lateral waves, but they fit each other at the NFFF boundary given by Eq. (9). Around 90°, 
the guide modes decay in the lower RI layers, but they are not covered in the discussion here because a guided mode 
is a bounded state that is equivalence to cylindrical wave and have a decay rate as l/i/Pi which means that its field 
is much larger than the results of LOA . 


4.2. The VED is in the middle layer 


(Az) 

Transmission Transmission 

z 

♦ VED 

d 


.. 

“1 

p 






Method A - Method B 


Figure 9: Physical meanings of SPs when a VED is located in the middle layer. Method A: blue lines; Method B: yellow lines. 
For transmissions in the top layer, two cases, both involving multi-reflection processes, are distinguished in Method A, where 
the propagation direction of the light is upward or downward after emitted from the VED. For a clear illustration, only the 
mode m=0 is demonstrated here. Method B sets the start point on the top interface right above the dipole. The transmissions 
in the bottom layer are similar to that in the top layer, but with the transmission direction reversed. 


Next, we discuss the case when the VED is located in the middle layer {d 2 < 0 < di), as shown in Fig. 9. For this 
configuration, the scattering coefficients Ci and Cq are expressed as 


Cl — ^21 


C 3 — t23 


gik2zdi 

f — (rfl — 

X — 7'2i7'23e^^^^^ 


r^^gik2z{di-2d2) 

X — 1*217‘23e^^^^^ Cl ”^ 2 ) ’ 
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X — 7*217‘23e^^^^^ Cl ~i^2) 
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For the \Ez\ in the top layer, Method A expands the two terms in Ci using geometric optics series, which is then 
substituted in Eq. (1) to get the SPs: 

{ 7 (m) 

p - - di) - -^[2m(di - d2) + di] = 0, 

p) ( 14 ) 

p - - di) - -^[2(m + l)(di - d2) - d2] = 0. 

^Izs ^2zs 

On the whole, Eq. (14) represents the light transmitting to the top layer after multi-reflections. The first equation 
represents the case of light propagating upward after being emitted from the VED. The second equation indicates the 
fact that light propagates downwards first, and after is reflected by the lower interface, it then propagates upwards. The 
intuitive picture is illustrated in Fig. 9. This course is also reflected by the expression of Ci in Eq. (13). Concerning 
the ray interpretation of SPs, we can write the asymptotic expressions similar to Eq. (6a). The transmission picture 
of Method B is the same as the above ones, and its asymptotic expressions are similar to Eq. (6b). The transmission 
to the bottom layer is similar to that to the top layer, but with the direction in reverse. 



0 (degree) 

Figure 10: \Ez\ as a function of angle 6 at distance r=100Ao for two configurations. The two critical angles are 9ci ~(180° — 
17°) = 163° and Od ~(180° — 26°)=154°. Structural schematics are illustrated in the insets. (|di|, |d 2 |, r) = (0.1,10,100)Ao. (a) 
[e : 1 — 2.25 — 12]. (b) [e : 2.25 — 12 — 1]. The gray area in (a) represents the second total internal reflection region. 


The \Ez\ distributions when the VED is in the middle layer of trilayer structures are shown in Eig. 10. Eigures 
10(a) and 7(b) have similar far-held distributions, and Figs. 10(b) and 8(b) as well. Therefore, the discussions there 
are valid for Fig. 10. The correctness of Eq. (9) is again verihed, where the parameter D is now the thickness of the 
middle layer. Moreover, Fig. 10 also shows that Method A provides more accurate results than Method B. 

It is known from the discussions above that conclusions about the NFFF boundary in a trilayer structure are similar 
to that obtained in a bilayer structure. The boundary in a higher RI layer is mainly determined by the lateral waves 
and satishes Eq. (9), while the boundary in the lowest RI layer is about ten wavelengths. 


5. GENERAL CONFIGURATIONS 

In this section, we generalize our conclusions to multilayer structures, and verify the universality of Eq. (9). 

Eor multilayer structures, light multi-reflects in each layer in the intermediate region. Consequently, the expressions 
of the field will have a recursive fashion, which will make the ray interpretations very complicated and limit the 
applicability of Method A. However, Method B is still of the simplicity in mathematics and is suitable for multilayered 
structures. 

The NFFF boundary in the multilayered structures is mainly affected by the lateral waves as well. The attenuation 
behaviors of these waves are similar to the ones shown in Fig. 4, so that Eq. (9) is correct in the order of magnitude. 
The value of the dimensional parameter D depends on whether the VED is in the top layer or in the intermediated 
region. When the VED is in the top layer, D equals the distance between the VED and the lowest interface at 
z = dq-i] when the VED is in one of the intermediate layer, D is the distance between the highest interface at z = di 
and lowest interface at z = dg-i, which is the total thickness of the intermediate region. The multi-reflection processes 
do not affect the D value for the following two reasons. Firstly, the multi-reflections have a quick convergence. All 
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the results of Method A shown in the figures converge when m is up to 2 ~ 4, which does not change the order of D. 
Secondly, as can be seen from Fig. 4, the lateral waves have a nearly r~^ decay when D is small, and have a more 
rapid decay when D is relatively larger. Therefore the influence from the multi-reflection may decay to a negligible 
one within the distance given by Eq. (9). Moreover, it is this decay characteristic of the lateral waves that makes the 
NFFF boundary almost independent of the wavelength. 

For other field components of the dipole radiations and the cases of different dipole orientations, we have numerically 
verified the correctness of Eq. (9). Besides, all these cases depend on the evaluations of the scattering coefficients 
related Sis, but the field expressions are slightly different, which results in different far-field patterns over the obser¬ 
vation angles. However, in the asymptotic analysis, they have the same SPs and similar asymptotic expressions and 
branch cut contributions. Thus, Eq. (9) is applicable. 

In the beginning of Sec. 2, we have assumed that all of the layers were lossless. Now let us discuss what about 
the case where the intermediate region is not lossless. In such a situation, two cases are distinguished to discuss the 
NFFF boundary: the loss is large or small. When the loss is large, as in undersea communications where short waves 
may not reach the sea bottom 22l-l24|. the air-ocean-earth trilayer model simplifies to an air-ocean bilayer one. From 
the discussions in Sec. 3, a distance of IOAq is enough to differentiate the near field and far field in air. When the 
loss is relatively small, the amplitude of light decreases as it arrives the bottom interface, which accordingly leads to 
a decrease of the differences between the asymptotic results and the accurate ones. Generally speaking, losses make 
the NFFF boundary reduced: its value will be less than that given by Eq. (9). 

Our conclusion can extend to bulk sources. On the one hand, within the scope of volume integral method Q, a bulk 
source can be considered as a superposition of dipoles. The far-field radiations of each dipole have been fully discussed 


the near fields of an arbitrary source 


above. On the other hand, within the scope of NFFF transformation (?. 27. 
can be converted to equivalence surface dipoles on a virtual closed surface. Therefore the dimensional parameter D 
should anchor to the brightest spots in the near field. Of course, our conclusions can be considered as an applicable 
scope of the NFFF transformation. 


6. CONCLUDING REMARKS 

In this paper, we have investigated the far-field asymptotic behaviors of dipole radiations in stratified backgrounds 
and obtain a universal empirical expression of NFFF boundary, Eq. (9). The asymptotic results are compared with 
the accurate numerical ones in various configurations to make sure the universality of Eq. (9). The NFFF boundary 
is mainly affected by the lateral waves, which correspond to branch point contributions in the Sis and are of a higher- 
order attenuation. In Eq. (9), the parameter D plays a key role, and it is much larger than the operating wavelength. 
As a result, the NFFF boundary in a stratified background is totally different from that in vacuum. To be more 
specific, Eq. (9) describes the boundary between the intermediate and far fields. However, since the intermediate 
field is vaguely defined in optics, we still call it as the NFFF boundary. In the case that the observation point is in a 
region where its RI is the lowest in the whole structure (usually air), the lateral wave decay rapidly, and the NFFF 
boundary is about ten wavelengths. 

It is believed that our conclusions are very helpful in understanding and applying the far-field approximation. In 
electromagnetic simulations, such as the finite-difference time-domain method and the finite element method, the 
far-field results are obtained by a NFFF transformation where the far field approximation, or the LOA, is employed. 
The NFFF boundary presented here reveals the applicability of the NFFF transformation, especially in the forbidden 
region. Moreover, we compare the different treatments for SPs in the asymptotic method and improve the accuracy 
according to the ray theory (Method A). The equivalence between the results of the asymptotic method and reciprocal 
theorem is demonstrated. 
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